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Abstract. We experimentally study the relaxation dynamics of a co- 
herently split one-dimensional Bose gas using matterwave interference. 
Measuring the full probability distributions of interference contrast re- 
veals the prethermalization of the system to a non-thermal steady state. 
To describe the evolution of noise and correlations we develop a semi- 
classical effective description that allows us to model the dynamics as 
a stochastic Ornstein-Uhlenbeck process. 



1 Introduction 

Non-equilibrium dynamics and relaxation processes in isolated quantum systems con- 
stitute a fundamental question in many fields of physics pQ . Open problems occur on 
vastly different energy, length and time scales, ranging from the dynamics in the early 
universe and heavy-ion collisions [2] to the subtle coherence and transport properties 
in condensed matter systems [3 4 5 6 7 8. However, in contrast to thermal equilibrium 
situations both the theoretical understanding of non-equilibrium phenomena as well 
as the experimental techniques to study complex transient states are still in their 
infancy. In particular, it remains an open question how an isolated quantum many- 
body system reaches thermal equilibrium from a given initial non-equilibrium state, a 
question touching the very fundamentals of quantum statistical mechanics |H9ll0llTj . 
This relaxation has recently been shown experimentally to proceed via an intermedi- 
ate non-thermal state |12j exhibiting prethermalization [5] , where the system already 
demonstrates certain features similar to thermal equilibrium, while it is still strongly 
non-equilibrium. In [12j , this decay was modelled using a Tomonaga Luttinger liq- 
uid formalism [13114] . Here we show that the correlations in the system during the 
relaxation can be described by an effective model based on an Ornstein-Uhlenbeck 
stochastic process [15] . 

The system we study is a repulsive, one-dimensional (ID) Bose gas in the quasi- 
condensate regime [16117] . In this regime density fluctuations are suppressed, and 
noise and correlations arc determined mainly by strong phase fluctuations along the 
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Fig. 1. Non-equilibrium dynamics of a coherently split ID Bose gas. An initial phase fluc- 
tuating f D Bose gas is rapidly split into two uncoupled gases, such that each gas has almost 
the same phase profile (j>i(z) ~ <f>2{z). Right after splitting, the relative phase profile 6{z) is 
thus almost flat, reflecting the strong correlations between the gases. During the evolution 
of the system these correlations are expected to be lost, leading to a randomization of the 
relative phase profile, which can be directly studied using the corresponding decrease of the 
longitudinally integrated interference contrast C(L) in time-of-flight. 



longitudinal axis of the system. This is in contrast to true three-dimensional Bose- 
Einstein condensates where long-range order allows the characterization of the order 
parameter by a single phase. As we will see in the following these phase fluctuations 
can lead to complex dynamics in ID Bose gases. 

To initialize the non-equilibrium dynamics, we rapidly split an initial single ID 
Bose gas coherently into two uncoupled gases with almost identical longitudinal phase 
profiles [13114118] . This situation is depicted in Fig. Jlk). While the two resulting gases 
contain all the thermal entropy of the initial gas in their common degrees of freedom 
(e.g. 9 com {z) = 4>i(z) + 4>2(z)), the relative phase profile 0(z) = 4>i(z) — 4*2( z ) is influ- 
enced only by very weak fluctuations resulting from the splitting process. Here, cf>i t 2 
denote the longitudinal phase profiles of the first and second gas after splitting, re- 
spectively. As these common and relative degrees of freedom are populated according 
to completely different energies, the system is expected to relax, eventually leading 
to a randomization of the relative phase profile 1191120121122) . The question we address 
here is how far this randomization will proceed. In other words, will the two gases 
completely forget their initial correlations as commonly associated with the approach 
to thermal equilibrium, or not? 

In our experiment the evolution can directly be observed by performing matter- 
wave interferometry in time-of-flight |23|24|25] , where the integrated contrast of the 
interference pattern is directly related to the relative in-situ phase profile as shown 
in Fig. ^p) . We will demonstrate in the following that this allows for a complete 
characterization of the transient states that are reached during the evolution. 



2 Experiment 

To prepare a single ID Bose gas, we use standard evaporative cooling in a magnetic 
microtrap on an atom chip The experimental setup is depicted in Fig. We typ- 
ically trap degenerate gases containing 2-10 x 10 3 Rb atoms in the F = 2,m,F = 2 
state. This corresponds to a peak linear density of rim ~ 20 — 80atoms/^m. The 
proximity to the current carrying structures of the atom chip allows the realization 
of very tight radial trapping such that the gas becomes very anisotropic and behaves 
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Fig. 2. Experimental setup, (a) A single ID Bose gas is realized on an atom chip. The 
atoms are initially trapped in a static magnetic trap formed about 100 /im below a micro- 
fabricated wire. Longitudinal confinement is provided by additional U-shaped wires, (b) 
The resultant harmonic trapping potential can be smoothly deformed into a double well 
trap with separation d ~ 3.5 /jm and barrier height Vb — 3 kHz by coupling the atomic 
Zeeman levels of the F — 2 manifold using radio-frequency radiation which is applied via 
parallel wires adjacent to the central trapping wire. This allows the realization of the initial 
non-equilibrium state described in Fig. ([T]). (c) Experimental chamber with the atom chip 
mounted upside down in the center. Two absorption imaging systems are used to probe the 
gas longitudinally (matterwave interference) and transversally (temperature measurement). 
The latter imaging system is also used for the optical pumping scheme described in Fig. pi). 



quasi-one-dimensional |17l30j . i.e. fcgT, fi < hui, where /i ~ 300— 1000 Hz is the chem- 
ical potential of the gas and ui ~ 2ir x 2 kHz is the typical radial trapping frequency. 
The temperature T ~ 20 — 150 nK is determined from measurements of the two-point 
density correlation function of the atoms obtained by imaging them transversally in 
time-of-fhght [27128129]. 

To realize the non-equilibrium experiment introduced above, we deform the static 
harmonic trap of the atoms into a fully controllable double well potential by applying 
radio-frequency radiation via additional wires on the atom chip [25131] . The typical 
timescale for the decoupling of the two gases during the fast splitting process is on the 
order of 0.5 ms. After the splitting, we hold the atoms for a variable evolution time 
in the double well potential. To study the resulting matterwave interference pattern 
between the two halves of the system, all trapping potentials are turned off and the two 
gases expand and overlap in time-of- flight. Information is extracted using absorption 
imaging along the longitudinal axis of the cloud |29j . This configuration automatically 
performs an integration of the interference pattern as depicted in Fig. ^p). To realize 
a certain integration length L in this process, we use an optical pumping scheme (see 
Fig. I3J) to spatially select atoms within a certain L. This procedure is technically 
simple to implement, very stable and avoids the complications connected to direct 
imaging of the interference pattern by reflection from the atom chip as discussed 
in [T2]. 

To study the dynamics we extract the interference contrast C{L) as a function 
of L and time. Typical interference pictures demonstrating the decay of contrast 
with time are shown in Fig. (|3ji) . This decay has been studied in detail in |12|25|38] , 
showing a rapid evolution over a timescale of approximately 10 ms, followed by the 
emergence of a steady state. This observation demonstrates that the relative phase 
does indeed partially randomize. However, it does not provide any information about 
the processes that are responsible for the decay or on the nature of the steady state. 

To extract more information about the dynamics we thus go beyond simple mean 
values and measure the full quantum mechanical probability distribution function 
(FDF) of the contrast C(L) |12|13)24|34"|35] . While the mean contrast can be con- 
nected to the two-point phase correlation function of the gas [16132133] , the FDFs 
contain also information about all higher even correlation functions and consequently 
allow the characterization of many-body quantum states with much more detail. Ex- 
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Fig. 3. Spatially selective matterwave interference, (a) To change the integration length L 
in the imaging process the cloud is transversally exposed to a 175 /ts- long optical pumping 
pulse after 1ms time-of- flight (TOF). This renders the edges of the cloud invisible to the 
imaging process. The desired length L is protected from the optical pumping by imaging the 
shadow of a variable sized target onto the center of the cloud, (b) Atomic level scheme for 
imaging and optical pumping. Imaging is performed on the F — 2 — > F' — 3 transition of the 
D2 line. For the optical pumping we employ the F = 2 — > F' — 1 transition, from where the 
atoms predominantly decay to the dark F = 1 state after scattering only very few photons, 
(c) Result of the optical pumping for various integration lengths L observed after 12 ms 
TOF expansion using the transversal absorption imaging system, (d) Typical longitudinal 
absorption pictures and corresponding line profiles for L = 110 /im and L = 22 fim after 
2.5 ms and 17.5 ms of evolution. The contrast C(L) is quantitatively extracted by fitting the 
line profiles with a periodically modulated Gaussian. The shortest L possible is determined 
by the remaining number of visible atoms, which decreases with decreasing L. 



perimentally, this has successfully also been used before to study ID gases in [24] 
and out of [12] thermal equilibrium, as well as the dynamics of an unstable quantum 
pendulum 36] . 

To extract the FDFs, we repeat the experiment many times with the same initial 
conditions and measure the different outcomes of the contrast. Note that this proce- 
dure is fundamentally different to the averaging performed in a 2D optical lattice |37] . 
There, typically many ID gases are realized and probed in parallel which means that 
only ensemble averages are accessible. Due to the central limit theorem the statistics 
of those average values is typically approximately Gaussian, meaning that the infor- 
mation contained in the higher moments (i.e. the special shape of the FDFs) is not 
accessible in this way. 

The results of our procedure are shown in Fig. Q . Initially we find peaked distri- 
butions on all length scales, directly reflecting the coherent initial state. The distri- 
butions develop markedly different forms for different L during the evolution. For an 
integration over the whole length of the cloud (L = 110 fim) the distributions quickly 
develop an exponential shape. This exponential shape is characteristic for a state 
where the correlation length is much shorter than L. On long length scales, our in- 
terferometer thus seems to have lost all its initial correlations. On the other hand, on 
short length scales (L = 22 /im) the distributions remain peaked for all times probed 
allthough L is still significantly larger than the thermal coherence length. This shows 
that some of the correlations of the initial state did indeed remain in the system. 

Similar to the mean contrast |12I38| , the FDFs change strongly during the first 
10 ms of the evolution, but then reach a steady state. We can use the FDFs to di- 
rectly show that the steady state the system reaches, is not thermal equilibrium. For 



Will be inserted by the editor 



(a) 



(b) 



110 |i m 




1 2 



60 jim 



1 2 



38 |xm 




27 jim 



22 jim 





>frL 






1 2 





2.5 ms 

T af| =14.1 nK 

J 8 =0.76 Hz 

erf 

4.5 ms 

T =14.1 nK 

eff 

J =0.24 Hz 

erf 

7.5 ms 

T =14.1 nK 
eft 

J =0.08 Hz 
eff 

12.5 ms 

T =14.1 nK 

eft 

J „=0.03 Hz 

eff 

17.5 ms 

T =14.1 nK 

eff 

J =0.02 Hz 
eff 

27.5 ms 

T =14.1 nK 

eff 

J =0.01 Hz 



equilibrium 

T=50 nK 



Fig. 4. (a) Full distribution functions (FDFs) of contrast C(L) 2 (histograms) revealing 
the underlying physics behind the dynamics of the non-equilibrium system. We find very 
good agreement with the effective Ohrnstein-Uhlenbeck model described in the text (solid 
line). Shaded areas denote the statistical uncertainty corresponding to the 170 samples 
used in the experimental data. The significantly different behavior on short and long length 
scales directly visualizes the multimode nature of ID Bose gases [38]. (b) Typical equilib- 
rium FDFs measured by splitting a thermal gas and creating two completely independent 
quasi-condensates. The strong difference in shape of the thermal FDFs and the observed 
non-equilibrium steady state directly demonstrate that this steady state is not thermal 
equilibrium. 



a system of two uncoupled quasi-condensates the thermal equilibrium state is char- 
acterized by the fact that both gases are completely independent. One can show that 
in this case relative and common degrees of freedom have identical temperatures, as 
expected for a thermal equilibrium state. We can prepare this situation in experiment 
by splitting a thermal, non-condensed gas using the double well potential. We then 
perform evaporative cooling in the double well, producing two completely independent 
quasi-condensates which do not have any knowledge of each other |24l33j . Mapping 
the distributions we find the characteristic exponential decay on all length scales 
probed (see Fig.[4b). This differs strongly from our observations for the steady-state 
of the non-equilibrium evolution of the coherently split system. We thus conclude, 
without the need to refer to any theoretical model that the observed state is not the 
true thermal equilibrium state of the system. 
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3 Understanding the dynamics 

We now introduce a model to describe the evolution of the system and the nature 
of the steady state, showing that the system reaches a prethermalized state. Pre- 
viously we have studied this system using a homogeneous, integrable Luttinger liq- 
uid description |12|13j , which has been shown to be a good approximation for the 
nearly-integrable dynamics of the trapped system. Here, we present a complemen- 
tary approach based on an effective stochastic model using an Ornstein-Uhlenbeck 
process 15 39 , which provides valuable insights into the correlation properties of 
the gases, is numerically very efficient, and allows us to take into account technical 
properties of the experimental setup more easily. This approach has been previously 
introduced to study equilibrium properties of ID 40J and 2D |4T] gases, providing a 
straightforward link between equilibrium and non-equilibrium dynamics in our ex- 
periment. 

A system of two uncoupled quasi-condensates can be described by the Hamiltonian 



H = I dz 



2 ' '^M^UW uu. 



E 



(i) 



where g — 2hoja is the ID interaction strength, m is the atomic mass, a ~ 5.23 nm is 
the 3D scattering length and \i = grim is the chemical potential. The index j = 1, 2 
denotes the two gases and nn> is the linear density in each one of t hem. We decompose 
the field operators into *Pj = y hj (z) exp[i(f>j (z)] . Here, 4>j(z) and hj(z) denote the 
operators describing the respective phase and linear density. This approach allows us 
to employ a generalization of standard Bogoliubov theory to describe the system, as 
density fluctuations in the two individual gases are suppressed and can be considered 
as small corrections [42143] . 

To describe the correlations between the two gases that are introduced in the 
splitting it is convenient to transfer the problem into symmetric and anti-symmetric 
coordinates which correspond to the common and relative degrees of freedom. We find 
that these degrees of freedom completely decouple in our system 1 14 . Consequently, to 
describe the evolution of the interference contrast, we only have to study the relative 
degrees of freedom, where 9(z) = <f>\{z) — <f>2(z) denotes the operator of the relative 
phase between the two quasi-condensates [18144) . The canonically conjugate variable 
is the relative number density v{z) = \fi\{z) — fi2(z)]/2. The expansion of 9(z) and 
i>{z) into plane waves reads 

*»=E^S' ^) = E^ £ 5- ( 2 ) 

fc/O V/L k^O V/L 

where we have assumed periodic boundary conditions with k = (2ir/£) x integer, C 
denoting the size of the system. Note that for the expansion coefficients [^L^fe'] = 

[O—kji'k'] — —iSkk'- In these variables the linearized Hamiltonian for the relative 
degrees of freedom is equivalent to a system of independent harmonic oscillators with 
momenta k and reads |14l32j 



£ = E 



(3) 



In the next step, we characterize the initial state of the dynamics. In the ex- 
periment the splitting is performed fast in comparison to the timescale set by the 
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inverse chemical potential £ sp ut < H/fi = £h/c, where £h = H/(mc) ~ 500 nm is the 
healing length and c = yj gum/m is the speed of sound. Thus there is no time for 
the atoms to correlate beyond £h along the length of the gases and the longitudi- 
nal density fluctuations are completely random. In particular, there is no correlation 
between modes with different momenta. Also, for each atom the decision of going 
to either one of the gases is random and uncorrelated, leading to a binomial distri- 
bution of atom number fluctuations in each small part of the ID system and the 
respective minimum-uncertainty relative phase distribution. Initially, at t = 0, the 
relative phase and density fluctuations are thus governed completely by the atomic 
shot noise and hence are given by (i>tu k )\ t —o = tiid/2 and (fi k 6 k )\f—Q — l/2fiirj. In 
contrast, the fluctuations of the ground state of the Hamiltonian (Eq. [3]) are given 

by (i>lA-)| G S = n lu S k /2 and (^ fc )| G s = l/2ni D 5fc, where S k = (1 + k 2 ^ 2 )~ 1/2 is 
the static structure factor. In other words, at t — we have a squeezed state: the 
relative phase fluctuations are suppressed with respect to their ground-state value. 
On the contrary, there is a huge excess noise of the number imbalance. 

In the subsequent evolution, this noise behaves semiclassically. In what follows, 
we substitute the operators v k and 9k with classical fields v k and 9k , respectively. The 
initial conditions for these classical fields are given by (|0 fe | 2 )| t=o = and (|^fc| 2 )|t=o = 
^id(1 — <5fc)/2. To obtain analytic results, we approximate the latter expression as 

/|„ |2v I _ / "ID/2, |fc| < f" 1 . , 

(M>|t=o-| Qj |fe| > e -i - ( 4 ) 

The Heisenberg equations of motion following from Eq. ([3| in the semiclassical limit 
are 

(hk 2 \ n • 1 fh 2 k 2 n 2 \ 

Vk = niD 7. — 9 k , 9 k = — h 2mc v k . (5) 

V 2m / ZwiiD V 2m / 



The solution yields 



„ , , h 2 k 2 /(2m) + 2mc 2 , . 

9 k {t) = ^(0 sinwfci, (6) 

n 1D nuJk 

where the eigenfrequency of a mode with momentum k is given by 
lfhk 2 \ fhk 2 mc 2 \ , / — — — 

Over time, the initially high density fluctuations periodically turn into phase fluctu- 
ations and dephase, which leads to the observed decrease in interference contrast in 
the course of the evolution. To describe the decay we calculate the two-point phase 
correlation function 

C( Z A = (Mz)ty(z)*Hz')Mz')) „ /eWM-w^h = e -m*)-°(z')f) ( 8 ) 
(\A(z)\ 2 )(\Mz>)\ 2 ) 

Writing t explicitly as an argument of C, Eq. Q can be expressed as 

rc0 dk 



C(z, t) = exp 







7T 



{\9 k {t)\ z ){l-coskz) 



(9) 



where z = z — z' and B k (t) is given by Eq. ^ . Recalling that significant semiclassical 
noise is present only in modes with \k\ < £T , we use the long- wavelength approxi- 
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mation for the mode eigenfrequencies: uik w ck. Then 
C(z, t) = exp 



, 5h „ sin 2 (ckt) ,„ , , 

2 / dk ^^-4(l-cosfcz) 



exp ■ 



2mc 2 t 
HJC 



z 



2ct 



S(\2ct - z\^) 



1 + 



2ct 



(10) 



where K, = 7rraiD£h is the Luttinger liquid parameter of the system, S(x) — (cosx — 
l)/x + Six and Six = J Q dyy~ x siny is the sine integral. 

For the sake of dcfinitcncss, we set z > from now on. If we are not interested in 
details of the correlation function behavior on the experimentally unresolved length 
scale < £hj we can simultaneously assume ct ^> 1, z^ 1 3> 1, and \2ct — z]^ 1 3> 1. 
In this limit we obtain 



C{z,t) 



exp[ 
exp[ 



- 7 rz/(2/Ca)] ! 
-irmc 2 t/(hlC)], 



0<z<2ct 
z > 2ct 



(11) 



We see from Eq. (Ill that the phase correlations first decrease exponentially and then 



stabilize at a certain level. The crossover between exponential decrease and constancy 
of C(z,t) and the value of C(oo,t) are time-dependent. 



The essence of our approach is now to to compare Eq. ( 11 ) to its counterpart C Gq (z) 



in the stationary regime of two tunnel-coupled quasi-condensates at finite temperature 
T(32] 



C cq (z) = exp 



-k t I. 



(l - e"^') 



(12) 



where At = h 2 niY>/mk^T = 1/kt is the thermal coherence length at temperature 
T characterizing the thermal fluctuations which lead to a random phase along the 
longitudinal axis. lj is the corresponding length scale of a coupling J which locks the 
relative phase of the two quasi-condensates [33] . This comparison leads to a conclusion 
that the time-dependent correlation function C(z,t) during the dephasing emulates 
that of a coupled equilibrium system with the effective temperature and effective 
time-dependent phase-locking length 



ZkB 



(13) 



A quasi-steady state is reached on a timescale t e 3> t, where the typical evolution 
time is [18] 

r = JC£h/nc = hJC/irmc 2 = nny^/c. (14) 

For the parameters of the dataset shown in Fig.Q we find r ~ 5 ms, in good 
agreement with our experimental observations that the steady state is reached for 
t e > 10 ms. The effective length scale is directly related to an effective coupling 
J e s(t) = /i/4mZj ff (i) 2 . In the limit t e ^> r the effective coupling vanishes and we 
recover exactly the prethermalized state predicted in [13] and observed in [12] . Al- 
though the system is still strongly non-equilibrium, it appears thermal-like and can 
be characterized by an effective temperature T e g, which is independent of the initial 
temperature of the gas. In particular, common and relative degrees of freedom are 
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indeed decoupled and have not yet equilibrated. All dynamics can be attributed to 
dephasing within the relative degrees of freedom p]3] . In the approach to the prether- 
malized state, lj plays the role of a characteristic length scale over which the system 
forgets the initial correlations for a given evolution time. The scaling of lj is that 
of a light cone |6l45j , reflecting that correlations in this many-body quantum system 
spread with a finite velocity given by the speed of sound. The typical correlation 
length reached in the prethermalized state is 

A eff = 2CT = 2/C&/7T = A T (T eff ). (15) 

For our typical parameters A c ff 3> Xt, showing again that the prethermalized state 
is characterized by a strong memory of the initial state, leading to much higher 
correlations than expected in thermal equilibrium [38] . 

Our effective model along with the assumption of Gaussian fluctuations enables us 
to simulate the FDFs for any evolution time t e and any sampling length greater then 
the healing length using an Ornstein-Uhlenbeck stochastic process [39] that develops 
in space, along the major axis of the trap, 

j- z <j>{z) = -r ] 1 ^{z) + f{z). (16) 

Here, f(z) is a random diffusion force characterized by (f(z)) = and (f(zi)f(z 2 )) = 
2kxS(zi — z 2 ), iy 1 describes friction. In our case, the ratio T e g/niD of the effective 
temperature to the linear atomic density defines the diffusion and the effective inter- 
well coupling provides a restoring force. The particular choice of correlation properties 



of f(z) is required to reproduce the correlation function given by Eq. ( 12 ) 39 40 and 



is in very good agreement with experiments |33j . The correlation function may differ 



from Eq. (12) on extremely short length scales < £h and hence f(z) deviates from 
white noise on this scale. However, this length scale is not optically resolvable in 
experiment and can therefore be neglected in the analysis. 

To construct the FDFs, we use an exact updating formula [46] to simulate in- 
dividual phase profiles. Properties of the experimental setup such as finite imaging 
resolution, trapping potential and time-of-flight expansion are directly taken into ac- 
count in this procedure |33) . A statistical analysis yields the FDFs. The method is 
found to be computationally fast and efficient, especially for finite-size atomic clouds. 

In Fig. Q we compare the results of the Ornstein-Uhlenbeck simulation to our 
experimental data. We find very good agreement for all length scales and evolution 
times probed, both for the non-equilibrium and equilibrium cases. 

All parameters have been independently extracted from experiment. Our approach 
directly reveals that the prethermalized state already shows thermal-like properties, 
while it is still strongly non-equilibrium. As predicted by our effective theory it can 
be described by an equilibrium distribution with an effective temperature. This tem- 
perature is determined only by the linear density and can be significantly lower than 
the true thermodynamical temperature describing the common degrees of freedom. In 
Fig. ([5]) we check the dependence of the effective temperature on the linear density of 
the gas, by fitting the experimentally obtained distributions with our model, leaving 
temperature as a free parameter. The results directly confirm the predicted linear 



scaling of Eq. ( 13 ), providing further evidence for our interpretation of the dynamics. 



4 Conclusion 

In conclusion, we have demonstrated prethermalization of a coherently split ID Bosc 
gas. While the near-integrability strongly suppresses true thermalization of the sys- 
tem, dephasing within the relative degrees of freedom leads to the establishment of 
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Fig. 5. Scaling properties of the effective temperature. The fast splitting process leads to a 
linear dependence of T e a on the density in each of the two gases, as given by Eq. (13 1. The 
datapoint at niu = 41 /im -1 corresponds to the dataset presented in Fig Q. The remaining 
datapoints have been adapted from [12] . 



an intermediate steady-state which has thermal-like features but is notably different 
from true thermal equilibrium. These dephasing dynamics can be well approximated 
by a semiclassical effective model based on an Ornstein-Uhlenbeck stochastic process. 
Our theoretical approach can straight forwardly be extended to the study of coupled 
gases, which could be used to reveal the dynamical properties of the Sine-Gordon 
model. Recent experimental advances now enable us to also image the full spatial 
structure of the interference pattern |12l38j , enabling, for example, the direct study 
of two-point phase correlations [33] . In the future, we plan to exploit this enhanced 
toolbox to shed light on outstanding questions about the second timescale of the evo- 
lution, where integrability is expected to be broken and thermal equilibrium is finally 
established. 
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